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Observations show that high-velocity jets stem from deeply embedded young 
stars, which may still be experiencing infall from their parent cloud cores. Yet 
theory predicts that, early in this buildup, any outgoing wind is trapped by 
incoming material of low angular momentum. As collapse continues and brings 
i/-) , in more rapidly rotating gas, the wind can eventually break out. Here we 

model this transition by following the motion of the shocked shell created by 
impact of the wind and a rotating, collapsing envelope. We first demonstrate, 
both analytically and numerically, that our previous, quasi-static solutions are 
dynamically unstable. Our present, fully time-dependent calculations include 
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cases both where the wind is driven back by infall to the stellar surface, and 
where it erupts as a true outflow. For the latter, we find that the time of 



breakout is 5 x 10 4 yr for wind speeds of 200 km s -1 . The reason for the delay is 
that the shocked material, including the swept-up infall, must be able to climb 
^ ! out of the star's gravitational potential well. 

We explore the critical wind speed necessary for breakout as a function of 
the mass transport rates in the wind and infall, as well as the cloud rotation rate 
Q and time since the start of infall. Breakout does occur for realistic parameter 
choices. The actual breakout times would change if we relaxed the assumption 
of perfect mixing between the wind and infall material. Our expanding shells 
do not exhibit the collimation of observed jets, but continue to expand laterally. 
To halt this expansion, the density in the envelope must fall off less steeply than 
in our model. 



Subject headings: circumstellar material — ISM: jets and outflows — stars: 
mass loss — stars: Pre-Main Sequence — Hydrodynamics: Shocks 
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1. Introduction 

The traditional picture of low-mass star formation in isolation (Shu, Adams, & Lizano 
1987) envisions an early period of pure accretion, prior to the appearance of any wind. This 
expectation is based on simple theoretical considerations. Since the star is building up mass 
from its parent cloud, the rate of mass infall exceeds that of any outflow. In addition, both 
theory and observation suggest that the speeds of infalling and wind gas are comparable. 
Thus, the ram presure from the collapsing envelope should crush the wind, preventing its 
escape from the stellar or inner disk surface. Why is it, then, that the most embedded stars, 
including those designated Class (Andre et al. 1993) are observed to produce vigorous 
winds? How do wind and infall simultaneously occur in a very young object? 

The answer, as has long been appreciated, is that neither flow is isotropic. Rotation 
and magnetic fields, in both the star and cloud, break spherical symmetry and alter the 
geometry of the wind and infall. Moreover, rotational distortion increases as the collapse 
of the parent cloud proceeds in an inside-out fashion (Shu 1977). Our goal in this paper is 
to examine numerically the interaction between wind and infall. Since the full situation is 
complex, it is best to proceed in stepwise fashion. We concentrate, therefore, on the effect 
of infall geometry, positing for simplicity a spherical wind. We also focus exclusively on 
the rotational influence. That is, we neglect any magnetic force on the infalling gas, under 
the assumption that decoupling from the ambient field has already occurred (Li & McKee 
2000). 

As in our first study of this topic (Wilkin & Stahler 1998; hereafter Paper I), we follow 
the dynamics of the thin shell formed by the colliding wind and infall. Paper I derived an 
evolutionary sequence of shells under the simplifying assumption of quasi-steady flow. That 
is, we took the expansion or contraction of the shell to be slow compared to either the wind 
or infall speeds. We now study the matter in more detail, solving the initial value problem 
of shell motion in a fully time-dependent manner. Our finding, in brief, is that the shell can 
either break out or recollapse to the star, depending on both the flow parameters and the 
epoch during cloud collapse. We derive the time of breakout as a function of the wind and 
infall parameters, and explain how to generalize this to the case of an anisotropic driving 
wind. 

In §2 below, we first detail the mathematical formulation of the problem, obtaining the 
time-dependent version of the dynamical equations. After describing our method of solution 
in §3, we show in §4 that our previous, quasi-steady shells are in fact dynamically unstable. 
We set up our initial conditions and nondimensionalization of the launch problem in §5, 
before presenting an overview of our numerical results in §6. We also rederive key findings 
in a heuristic, analytical fashion. Finally, Section 7 discusses the broader astrophysical 



- 3- 



implications of our study. In particular, an observational signature of the trapped wind 
phase may be fluctuations of the accretion luminosity, as an oscillatory shock structure is 
expected to develop prior to breakout of the wind. 

2. Formulation of the Problem 

2.1. Thin-Shell Approximation 

The supersonic collision of wind and infall leads to the formation of both an inner and 
outer shock front. At the wind speeds appropriate for low- mass stars (typically less than 
300 km s^ 1 ), the cooling behind these shocks is relatively efficient. (Recall the discussion 
in §2.2 of Paper I.) Hence we may describe the intershock region as a cold, thin shell, 
characterized by its radius R{6,i) and mass surface density, a(9,t). Here, 9 is the polar 
angle measured from the rotation axis of the parent cloud. As in Paper I, we assume this 
shell to be axisymmetric, i.e., invariant with respect to the azimuthal angle 0. 

Within the shell, two fluids with very different temperature and velocity come into 
contact. The internal shearing layers are subject to the Kelvin-Helmholtz instability, which 
leads quickly to turbulent mixing. We shall assume that the mixing is so efficient that we 
may describe the shell as a single fluid with a time-averaged velocity at each point (R, 9, t). 
This velocity has an azimuthal component arising from rotation of the infalling envelope. 
It also has a meridional component along the shell. Furthermore, we do not neglect, as in 
Paper I, the velocity of the shell normal to the surface. 

The evolution of the shell is governed by the incident fluxes of mass and momentum 
from the wind and infall. It is also influenced critically by the gravitational attraction of 
the star for the shell material. The central difference of the present study and Paper I is 
that we no longer assume the evolution to be quasi-steady. That is, the properties of the 
shell are allowed to change over the time required for material to travel along the surface. 

2.2. Description of Infall and Outflow 

The infall arises from the gravitational collapse of a dense cloud core within a larger 
molecular cloud. To obtain the resulting accretion flow, we idealize this core as a singular 
isothermal sphere rotating rigidly with angular velocity Q - Our calculation does not 
address the formation of the core itself, and neglects any turbulence in the initial structure. 
Recent theoretical work (e.g., Cho, Lazarian, & Vishniac 2002; Vazquez-Semadeni, 
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Ballesteros-Paredes, & Klessen 2003) suggests how cores might condense out of a turbulent 
flow. Nevertheless, observations continue to show that nonthermal motion in the core itself 
is relatively small (Barranco & Goodman 1998). Once the infalling gas becomes supersonic, 
the effects of this turbulent component are expected to be minor. Our adoption of a 
spherical body, while clearly an idealization, should be of little consequence for the flow 
deep within the core's center. 

In any event, the distributions of infalling density pi and velocity u, depend on time, 
because the collapsing region spreads out as a rarefaction wave at the isothermal sound 
speed a , gradually engulfing material of higher specific angular momentum (Shu 1977). We 
adopt the infall model of Cassen & Moosman (1981), which represents the inner limit of 
the full collapse. This limit applies to radii that are small compared to the expansion front 
Rexp = flo t, where t is the time since the collapse began. 

We may gauge the total rate of infall by examining the transport of mass across a 
surface located well inside the expansion front. According to the inside-out collapse model 
of Shu (1977), this transport rate is: 



where m = 0.975. The infall itself is characterized by a time-dependent length scale, the 
centrifugal radius R ce n- This is given by 



Our reference value for Q is chosen to be consistent with inferred rotation rates based 
upon velocity gradients in cloud cores lacking embedded IRAS sources (Jijina, Myers, & 
Adams 1999). For an equator-on rotating core, this corresponds to a velocity gradient of 
0.6 km s -1 pc _1 . In regions where r R cen , rotational distortion of the flow is significant, 
and some of the infalling matter strikes the equatorial plane (8 = vr/2), forming a 
circumstellar disk. When r » R cen , the infall is nearly spherically symmetric. Equation 
(1) gives the sum of the mass per unit time impacting the star directly and that entering 
the disk. 

In terms of the nondimensional radial variable ( = R cen /r, the infall velocity 
components and density are (Terebey, Shu & Cassen 1984) 






(3) 
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,'GM*\ 1/2 /cos# D -cos#\ / cos#\ 1/2 
u i,e = — —~a 1 + -Z^TT (4) 



Uirh = 




/ GM* \ 1 / 2 sin 9 Q 

^inTl 1 l5j 



= - i ^[l + 2CP 2 (oos6l )r 1 . (6) 

The function P 2 (cos^ ) is the Legendre polynomial, where # is a Lagrangian variable that 
labels the streamlines. Specifically, 9 Q is the initial polar angle of each fluid element, just 
before it is overtaken by the rarefaction wave. This angle is given implicitly in terms of the 
instantaneous coordinates ( and 9 by the trajectory equation 

cos#o — cos# 

C= sin 2 9 cos 9 ' ( ' 

Our goal is to assess the collimating influence of the anisotropic infall. Accordingly, we 
idealize the wind to be spherically symmetric, with an associated mass transport rate of 
M w . If the star lies at the origin of a spherical coordinate system, then the wind density at 
a radial distance r is 

M w 

Pw = A^V W (8) 
Here, V w = u WyT is the wind velocity. Note that the components u Wt e and u w ^ are both zero. 

We do not investigate the structure of the star, which is taken simply to be a spherical, 
gravitating mass. The value of this mass is 

M* = Mit. (9) 

Thus, we assume that both the matter entering the disk and that colliding with the infall 
are efficiently recycled to the star. How gas within the shocked shell reaches the disk is, 
of course, one of our main concerns. On the other hand, the transport of matter within 
the disk through internal torques is beyond the scope of this project. Once the shell truly 
moves dynamically, mass accretion onto the star and disk is halted, and we freeze the stellar 
mass at the value corresponding to the launch time of the shell (See §4). 



2.3. Time-Dependent Equations: Mathematical Derivation 



As in Paper I, we consider a small, three-dimensional patch of the shell, whose center 
is located at (R, 9, 0) within a global, spherical coordinate system centered on the star. We 
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also utilize a global, Cartesian system (x,y,z). (See Figure 3 of Paper I, i.e., Figure 1.3) 
To help in following our necessarily abbreviated derivation of the equations, we advise the 
reader to consult §2.3 of Paper I. 

Figure 1, a slightly modified version of Figure 1.4, shows the representative patch in 
more detail. We let 7 denote the surface tilt, i.e., the angle between the patch normal n 
and the radial direction from the star r. The upper and lower faces, i.e., those depicted 
with the largest areas, coincide with the inner and outer shock fronts. The narrow faces to 
the left and right of the upper one trace loci of constant polar angle 9, while the remaining 
two have fixed <fi. Note that the length As is a small increment of the global coordinate 
s measuring distance along the shell from the pole to the equatorial plane. Inspection of 
Figure 1 shows that As = Rsecy A 9 and that the patch width is Aw = Rsin9 A0. The 
surface area of either the upper or lower face is AsAw = AA9A(f>, and the total patch 
volume is AsAwAn = AA9A<pAn. Here An (called Ah in Paper I) is the patch thickness. 
The factor A is given by 

A = J R 2 sin0sec7. (10) 

For a mathematical description of the shell evolution, we let R = R(9, <j), t) be a 
dependent variable. Then the tilt angle 7 is given by tan 7 = —R'/R, where R' = dR/89. 
During the evolution, we want the sides of our patch to retain fixed angular positions 
in 9 and 0. On the other hand, we allow the upper and lower surfaces to move in and 
out radially, in order to follow expansion or contraction of the shell. Let us denote by 
(u r , Ue, u^) the velocity of matter within the shell. (Note the lack of additional subscripts, 
which we use to denote wind or infall.) The velocity component of this fluid pointing along 
the patch normal is 

u n = M r cos7 + sin 7. (11) 
Now shell expansion gives the patch itself a normal velocity 

v n = -^cos'j. (12) 

Our kinematical constraint on the patch motion is simply u n = v n . Combining equations 
(11) and (12), we find 

dR 

— = Ur J rUd tan7. (13) 

Proceeding to the dynamical equations, we first derive an expression for AQ, the rate 
of change of any physical quantity within the patch. Let q denote the volume density of 
this quantity. This density may change with time, as may the patch volume. We thus have 

AQ = ^(qAnAsAw^J = (AqAn^J A0A0. (14) 
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where the partial derivative is at fixed 6 and 0. 

The change in our quantity comes in part from advection into and out of the patch. 
The rate here depends on both the wind and infall velocities, and those within the patch 
itself. Note that the relevant wind and infall velocities are the normal components relative 
to the patch, since the latter has the normal velocity u n . Let q w and g« be the density 
of the quantity in the wind and infall, respectively. Then these flows make an external 
contribution to AQ: 

AQ cxt = [iu w . n - u n ) q w - {u it n - u n ) qi] As Aw, 

= ^K,n?» _ U 'i,n(li} A9A(f) . (15) 

where we have used primes for the relative normal velocities. 

We must also account the internal contribution to advection, i.e., the effect of flow 
within the shell. Part of this flow is due to motion in the ^-direction. Noting that the radial 
thickness of the shell is Ar = (An) sec 7, Figure 1 shows that 

AQ inte = [u e qAr Aw]e-&o/2 - [u e q Ar Aw] e+Ae/2 
The final contribution to advection stems from the 0-velocity Here we find 



AQ int ,s = [u^qAsA n]^ M/2 - [« q A s A n}^ +A ^ /2 

dq 

= —us R sec 7 — A n A 6 A 0. 

0(J) 

Apart from notation, this equation is identical to equation (1.18). 



(17) 



We now sum all the advective terms to form AQ adv . Combining equations (15)-(17), 
we have 



A Qadv 



d_ 
"ffl 



A^qAn ) -+ A 



U w,n1w U i,n ( li 



us d 

— qt(? An) 
w o<p 



A9A< 



(18) 



where w = R sin 9 is the cylindrical radius. 



In order to express mass conservation, we let q be the mass density p. Note that this 
quantity, like the others we shall be considering, is formally infinite in the limit An — > 0, 
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whereas p An is the finite surface density a. Since there are no sources or sinks of mass, we 
demand that AQ = AQ adv . We find 

d ( . \ d ( .ouq s 



dt 



Aa ) + w( A ^jf) =A {P*> U ",«-Pi<,n}- (19) 



Turning to momentum, we first note that the magnitude of the gravitational force on 
the patch is AF g = GM*aAA9A(f)/R 2 . This force would increase the momentum within 
the patch even if there were no advection. Suppose we choose our patch to be located at 
4> = in the global, spherical coordinate system. Then the Cartesian components of the 
gravitational force become AF 9iX = —AF g sin9, AF giV = 0, and AF 9yZ = —AF g cos9. The 
physical conservation law is most directly expressed in terms of Cartesian coordinates. In 
the x-momentum equation, we let q = pu x in equations (14) and (18), and demand that 
AQ = AQ adv + AF gx . For the patch centered on = 0, u x ~ u m and du x /d(fi ~ — u^. 
Here u m is the cylindrical radial component of velocity. We obtain 

d ( . \ , d ( .ou e u m 



dt 



, (T a,, GM,a . 
Pi u itn Ui tW + ^— sin^. (20) 



The y-component of momentum conservation is simpler since AF M = 0. We now let 

q = pu y , and use u y ~ and du y / d<p ~ u m to find 

d ( . \ d ( .au e us\ f . 

oi [ A ° u * ) + oe [ A -r- J =A \ pw Vn <W 

-PiU t Uij \. (21) 

W ) 

We may obtain a simpler form of this 0-force equation by instead writing it in terms of 
the z-component of angular momentum. First note that using equation (13), we have the 
relation dw /dt + {uq / R)dzu / '89 = Multiplying this equation by Aau^, and adding it 
to w times equation (21), we obtain 



-PiU'i,n U i,4>)- ( 22 ) 



Similarly, the z-force equation is obtained using q = u z 



d ( . \ d ( a u e u z , . . 
— [Aau z j + — \A— — ) =A\ Pw u w ^ n u W)Z 



GM*a i 
-p i u itn u itZ ^—cos9\. (23) 
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We have seen that beginning with Cartesian components of velocity, the azimuthal 
symmetry of the problem leads naturally to equations in terms of the cylindrical components 
of velocity. However, for a radial driving wind, it will be most useful to use spherical 
components. By taking linear combinations of equations (20) and (23), and substituting 
u-oj = u r sin 9 + ue cos 9 and u z = u r cos 9 — uq sin 9, we may recast them in terms of the 
velocities u T and u e : 

-p iU ' itnUitr +-Jt — - , (24) 



R R 2 
d ( A \ d ( A <7UqU0 

dt 



A aUe ) + 1)6 { A = A {p™ U 'w,n u w ,e 



auicotO] 

-PiU i:n Ui, e + — ^ 1. (25) 

Previously, Giuliani (1982) derived an equivalent set of equations, including magnetic fields, 
but without accounting for gravity or rotation. 



2.4. Time-Dependent Equations: Lagrangian Form 

In practice, it is simplest to solve the evolutionary equations by recasting them in 
Lagrangian form. Consider first the comoving derivative of any quantity Q within the shell. 
In the Eulerian description we have used until now, Q is a function of t, 6, and <fi. It also 
depends implicitly on radius, since, at fixed 9 and <f>, the radius R varies with time. The 
comoving derivative is therefore 

DQ _ dQ u e 0Q dQ 

Dt ~ dt R 89 R sin# d<p ' 1 ' 

We now let Q be, in turn, the radial and angular positions of a moving fluid element that 
constitutes part of the shell. This element is no longer constrained in 9 and 0, as was our 
patch. Retaining the old notation for these coordinates, now serving as dependent variables, 
we substitute into equation (26) to find 



DR 

~Dt ~ Mr ' 
D9 u e 

Ut ~ ~R' 
D<P u<t> 

~~Dt ~ Rsin9' 



(27) 
(28) 
(29) 
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In deriving equation (27), we have used both equation (13) and the definition of 7 in terms 
of R'. Equations (27)- (29) describe both the expansion of the shell and tangential motion 
along its surface. We next substitute for Q in equation (26) the quantity Ao. After utilizing 
equation (19), we arrive at the Lagrangian expression for mass conservation, which we may 
write as 

DIMAo)) , , d (u e \ 

jy t = [p w u Wtn -p l u i J/(T-—^—j. (30) 

To obtain the r- and 9- components of the force equation, as well as the expression for 
angular momentum conservation, we substitute for Q the quantities u r , ug, and zuu^, 
respectively. We then use our Eulerian conservation laws and equation (30) in each case to 
find 

Du r ug I 2 sin 9 G M* 

~JJ~j~~ = lP™ U w,n U w,r ~ PiUi, n U iir \/a+ — + — g — , (31) 

Dug , , / / 1/ U r Ug I 2 COS 9 

= [Pw u wfi - piU itn u it g\/a — + ^ 3 , (32) 

jfi = [Pw u ' w ,Jw- Pi u 'i,J[\l °- ( 33 ) 

Here we have let I = zuu^ be the z-component of specific angular momentum, and have 
further defined 

C = ro K,* _lt ^ (34) 
l\ = w(ui^-u^). (35) 



Equations (27)-(29) and (31)-(33) may be summarized in vector form as 
DR 



Dt 



(36) 



D u . . , ... GM* A 

= [PwU^ n u w -piU^u^/a- -£2-r. (37) 

As expected, the acceleration of the fluid element depends on both the radial force of 
gravity and on the input of momentum from wind and infall. Note finally that the mass 
conservation equation (30) contains an Eulerian derivative on the righthand side, so that 
our equations are not in purely Lagrangian form. Indeed, the normal components of wind 
and infall speed that enter equations (30)- (33) also require, through the angle 7, the 
Eulerian derivative dR/d6. In practice, this mixed character of the equations, and the fact 
that 6 acts as both a dependent and independent variable, present no special difficulty for 
integration. 
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3. Method of Solution 

Equations (27)-(33) fully describe the motion of a fluid element within the shell. 
Because the latter is axisymmetric, there is no need to track the coordinate, and we may 
ignore equation (29). We treat the others effectively as ordinary differential equations in 
time, and use them to follow a collection of 50 discrete points that represents our shell. We 
evaluate the cross derivatives dR/d6 and d(u e /R)/d9, by numerically differencing R and 
ug/R. This method has been employed previously for the similar problem of steady-state, 
non-axisymmetric bow shocks and expanding supershells (Mac Low and McCray 1988; 
Bandiera 1993), and works well provided the spacing between adjacent trajectories is fine 
enough to determine numerical derivatives. 

For each of the 50 points, we integrated the equations in time using a fifth-order 
Runge-Kutta scheme. Cross derivatives were obtained at each point by using the two 
nearest neighbors, and were thus second-order accurate. Since the points move in a manner 
dictated by the local fluid velocity, they are unevenly spaced in angle. In practice, we 
approximated the shell segment by fitting a circular arc through each triad of points. The 
cross derivatives at the middle point were then obtained analytically from this curve. 

As the shell evolved in time, some of its representative points inevitably drifted toward 
the equatorial plane. In addition, points could move so close together that there was little 
variation in the intervening space. In either circumstance, we removed a point from the 
original set. We immediately replaced it with another, which we introduced at the part of 
the shell that was most sparsely covered. We always retained one point on the symmetry 
axis, where the boundary conditions R' = u e = = were applied. Here the mass 
equation was written in modified form to account for the vanishing of the sin 9 factor within 
A. 

We tested our code on three problems with known, analytic solutions. They include: 
(1) an expanding, spherical shell driven by an isotropic wind, (2) a shell driven by an 
angle-dependent wind within an r~ 2 density distribution (Shu et al. 1991), and (3) the 
bowshock created by the wind from a star that is moving into a uniform medium. The 
first two tests confirmed the code to better than 1% accuracy in a non-steady situation. 
For the last problem, Wilkin (1996) found analytic solutions for the steady-state shell 
configuration, including the surface density and flow velocity along the shell. Figure 2 
shows the expansion of our time-dependent bowshock as it approaches the steady-state 
endpoint. We have checked that the approach to equilibrium agrees with the calculations 
of Giuliani (1982). Figure 3 compares a number of quantities in the analytic solution with 
our numerical results, at a time when equilibrium has been reached. The fractional errors 
are between 10~ 2 and 10~ 4 , except close to the pole. Judging from the three tests, we 



believe our code is capable of tracking all variables in the time-dependent shell to within an 
accuracy of 1-2 percent. 



4. Instability of the Quasi-Steady Solution 



Before displaying our fully time-dependent results, we first discuss the quasi-steady 
solutions obtained in Paper I. 1 These shells are in dynamical equilibrium, with gravity and 
the infall ram pressure opposing the outward ram pressure from the wind. Because gravity 
plays a dominant role, we suspected that the shells might be dynamically unstable. (See §4 
of Paper I.) We now demonstrate this fact both analytically and numerically. 



We first examine the stability in the region of the polar axis. This approach extracts the 
basic physics while bypassing a full modal analysis, such as that done for wind bowshocks 
by Dgani, Van Buren, & Noriega- Crespo (1996). The total radial force per unit area acting 
on the shell near the axis is 



Note that we have omitted centrifugal terms, which are vanishingly small near the pole. 

The condition of normal force balance is that F r (R Q ) = 0, where R Q is the shell's 
equilibrium polar radius. We now imagine perturbing the shell, and examine the resulting 
change in F r . For this purpose, we consider the lowest-order oscillation mode of the shell, 
i.e., its breathing mode. We find the altered wind density p w from equation (8), and the 
infall density pi and velocity u i)T from equations (6) and (3), respectively. For the latter 
two quantities, we specialize to the pole by setting 9 a = = 0. To obtain the perturbed 
value of the mass density a, we ignore any additional mass swept up by the shell during 
its oscillation. Thus, to conserve mass, we write a = a Rl/R 2 , where a Q is the equilibrium 
value. Equation (38) then becomes 



4.1. Analytical Argument 



GM*a 



(38) 



F r 



M W V W 
AuR 2 



AttR 2 \ R ) 1 + 2C 



Mi {2GM*\ 1 / 2 1 



GM* a Rl 
R 4 



(39) 



: We draw the reader's attention to two typographical errors in Paper I. Equation (70) should have an 
overall minus sign on the right hand side, and equation (72) should read e = a sin 9 — r $^ l /i?$ t cos 7 . 
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where ( = R ccn /R. Note that we have also neglected alterations to the surface density due 
to the perturbation in the tangential motion. We showed in Paper I that such motion is 
generally small compared to the wind and infall speeds. 

We next nondimensionalize equation (39) by dividing through by M w V w /AitR 2 . Again 
employing subscripts for equilibrium values, we have 



R 



-2 



~fi 



R\-5/2 1 + 2 ( 



- f 

R J \R J 1 + 2 C V R 



R 



(40) 



The constants fi and f g are the nondimensional (fractional) contributions of the infall and 
gravitational forces, respectively, at the equilibrium position. These are given by 



fi 



h = 



Mi /2GMA 1/2 1 



M w \V*Rj l + 2Co' 
4 7T GM* a 

M W V W 



(41) 
(42) 



The requirement that R be an equilibrium radius means that f g — 1 — fi. Thus, we 
may eliminate f g from equation (40). Writing R/ R — 1 + 5 and recalling that ( oc l/R, we 
linearize in 5 to find the small radial force due to the oscillation: 



2-h 



3 + lOCo 



2(1 + 2 Co) 



5. 



(43) 



We see that there is a critical value for 



fi,crit 



4(1 + 2 Co) 
3 + lOCo 



(44) 



If the actual fi in the equilibrium solution exceeds fi tCr it, then the radial force is directed 
opposite to the displacement, i.e., the shell is stable (or overstable). Conversely, fi < fi iCr u 
implies dynamic instability. The minimum value of fi tCr u, that for infinite Co, is 0.8. That 
is, at least 80 percent of the confining force must come from the infall ram pressure, rather 
than gravity, to ensure stability. 



We now return to the exact results of Paper I. There we found fi to be 

fi = 



l | 3 (l + 2Co)(l + a(l + 2Co)) 2 

4 Co + \fQ 



3C/8 



(45) 



where a = M w /Mi is the ratio of wind to infall mass transfer rates. This result is taken 
from equation (1.60), where the three terms of that equation correspond, from left to right, 
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to the normal force due to wind, infall, and gravity, respectively. For a fixed value of Co, 
decreasing a increases /j. Thus, we consider the limit a — * 0, to obtain the largest possible 
value of /j. If we set a = in equation (45), then /j monotonically increases with ( . The 
maximum of fi occurs in the limit ( — > oo, a — > 0, in which case /j = 4/7. Since this is 
substantially less than the minimum value of f\ necessary for stability, we conclude that 
there are no combinations of a and ( a such that the shell is stable. While this analysis 
is based upon conditions near the symmetry axis, the small contribution of centrifugal 
effects at larger angles (Section 1.3.3) suggests that the same result should hold throughout 
the shell. One might argue that the steady-state solution might nevertheless be achieved 
in some average sense. While this has been found to be true for isothermal bowshocks 
(Blondin & Koerwer 1998; Raga et al. 1997), the calculations we next describe indicate that 
this is not the case in our problem. 



4.2. Numerical Demonstration 

We next use our time- dependent code to evolve shells starting near the steady-state 
conditions derived in Paper I. The latter calculation gives us the run of R, a, u t , for a grid 
of values. We first verified the validity of the steady-state solutions, checking that the 
net normal force on a fluid element is equal to the normal component of the centrifugal force 
it experiences. Despite this balance, the shells do not stay in their initial configuration, but 
evolve quickly. Shells begun at a size slightly larger than equilibrium expand, while those 
of slightly smaller size collapse, as expected by the analytical argument. A representative 
example is shown in Figure 4. The innermost two curves are polar radii for shells distorted 
slightly inward from equilibrium, while the outermost two curves represent shells begun 
slightly beyond the equilibrium radius. All other quantities in the shell, such as the run of 
surface density with polar angle, were initially the same as for steady-state. The equilibrium 
solution shown is that of the critical (innermost) model for a = 0.1 shown in Figure 1.9. The 
slightly shrunken, collapsing shell quickly accelerates to free-fall speed, while the expanding 
structure eventually decelerates at large radius as a momentum-conserving snowplow. The 
analytic argument of the preceeding section, together with these and similar numerical 
experiments, show convincingly that the quasi-steady shells are dynamically unstable. 
This basic result motivates our current study, where we consider the shell evolution as an 
explicitly time-dependent, initial value problem. 
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5. Time-Dependent Evolution: Basic Considerations 

5.1. Initial Conditions 

We now wish to follow the motion of the shell from its launch near the protostellar 
surface. Our shell is initially spherical and massless, with radius equal to that of the 
protostar, R*. As the structure evolves, it receives mass from both the wind and infall. 
Thus, even at launch, there are well-defined values for all the physical quantities of interest. 
Note, however, that the momentum equation (37) is singular if a vanishes. The initial 
velocity is found by setting the corresponding numerator to zero: 

Pw (u™ - Uo) - pi u' i n (Ui - Uo) = 0. (46) 

Here we have explicitly written the vector velocity differences. We define A = —piu' i n / p w u' w n , 
which is the ratio of mass fluxes from the wind and infall onto the shell. Then equation 
(46) becomes 

Uo = ^±^i. (47) 
1 + A 

To find A explicitly, we take the normal component of equation (46), obtaining 
Pw u w,n = Pi u ?,n- Defining rj = Pi/p w and taking the square root of both sides, we have 

(U w , n ~ Uo,n) = 1 1/2 (V - U i>n ). (48) 

In obtaining equation (48), we have accounted for the fact that u itTl < u Qjn < u w>n . Equation 
(48) implies 

A = (49) 

which yields the final result 

U ° = 1+^ ' (50) 

In equation (50), both Uj and r\ depend on 9. The normal component of this equation 
represents the balance of wind and infall ram pressures in the frame of the shell. This 
result is a generalization of the well-known formula for the speed of the planar bow shock 
driven by a steady jet (e.g., Blandford, Begelman & Rees 1984). In our case, however, the 
"ambient" matter is in nonuniform motion with velocity Uj. For a dense wind r\ <C 1, the 
shell tends to the wind speed, while for a rarefied wind with rj 3> 1, the shell collapses at 
the infall speed. 

To derive the initial rate at which mass is being swept up, we use equation (30): 

Da da . 
Dt ot 

= PwU' w , n - PiU'i,n, ( 51 ) 
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where we have dropped terms from equation (30) that vanish with a. The righthand side of 
equation (51) is to be evaluated using equation (50) for the shell's initial velocity. 

Our calculation requires values for the stellar radius i?*, the wind speed V w = \u w \, 
and mass loss rate M w . The distribution of infalling matter is specified by t, the time since 
the start of collapse, together with the sound speed a a and rotation rate fl Q of the parent 
dense core. Thus, there are six dimensional parameters. In principle, the protostellar radius 
is not independent of the others, but should be obtained from stellar evolution theory (e.g., 
Stahler 1988). However, we will simply adopt a characteristic protostellar radius and show 
how the results scale with this and with other choices of the dimensional parameters. 



5.2. Nondimensional Parameters and Equations 

To reduce the number of runs needed to explore parameter space, we cast the equations 
in nondimensional form. There are two time scales of interest - an evolutionary and a 
dynamical one. The first is the time over which the infall itself changes appreciably. In 
terms of the dimensional parameters listed previously, we define 

^fag) ■ (52 > 

According to equation (10) of Stahler et al. (1994), t ev is, to within a factor of order unity, 
the time after start of collapse when the centrifugal radius associated with infall attains the 
value i?*. 

We find the dynamical time by first assigning units of length and velocity. The first is 
i?*, and the second is the Keplerian speed at the stellar surface. Since the characteristic 
stellar mass is Mit ev , we have 

a 4 / 3 

V * = Witf*' (53) 

where we have left out the factor of m entering equation (1) for Mj. The dynamical time 
is now defined as or 

r.4/3 q i/3 

^dyn — 4/3 • \<J^J 

Oo 

Note finally that the mass of the shell, after a dynamical time, is of magnitude Mit dyn . We 
formally define our mass unit as 



a 



3 



R* OA 4/3 



M *» = 4^li^) ■ < 55 > 
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Here again we have left out the factor of m , and have introduced a factor of Air for later 
convenience. 

All physical variables can be written as dimensionless factors times appropriate 
combinations of i?*, td yn , and M dyn - Additionally, certain nondimensional ratios enter the 
final equations. One involves the mass transport rates in the wind and infall: 

- M w 

a = . 56 

Mi 

In this study, we only consider values of a less than unity. A second ratio is the wind speed 
divided by our dynamical unit of velocity: 

(57) 

Assigning canonical values of the dimensional parameters a a = 0.2 km s -1 , fi D = 2 x 10 -14 s -1 , 
and R* = 3i? , and a reference wind speed of 200 km s -1 , the reader may verify that v is 
of order unity. The third key nondimensional ratio is the initial launch time of the shell in 
terms of t ev : 

To = ti n it/t ev . (58) 

For r D < 1, the centrifugal radius has not yet reached the protostellar radius, and the infall 
conditions are those of nearly spherically-symmetric free-fall. For r D ^> 1, the initial launch 
is within the inner, anisotropic region of the infall. In order to compute the evolution of the 
shocked shell in dimensionless units, only the three nondimensional ratios a, is, and r need 
to be specified. 

Another obvious ratio is that of the dynamical and evolutionary timescales: 

tdyn 



+ 1 



(59) 



This quantity is of order 10~ 7 . Since our calculations span only dynamical times, e does not 
appear explicitly in the final, nondimensional equations (see below). Note that our previous 
physical units can be written as 

V, = e-^ao, (60a) 
t dyn = e^n~\ (60b) 
M dyn = e V3_!_ (60c) 
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The nondimensional equations include various quantities related to the infall. These, 
in turn, depend on the time t since the start of collapse. Since our basic temporal unit is 
tdyn, we may write 

t = e~ 1 Tt dyn . (61) 

Here, r is the nondimensional evolutionary time, written in terms of t ev . In practice, the 
difference between r and t is exceedingly small. Even if we were to evolve the shell for 10 4 
dynamical times, r would only change fractionally by order 1CT 3 . Thus, in the dynamical 
equations we treat r as a constant, and drop the unnecessary subscript on r D . 

To complete our description of the nondimensional problem, we give the dimensionless 
forms of the remaining infall and wind quantities. We write the infall vector velocity in 
terms of the Keplerian speed as Uj = uk f , where uk = \JG M*/R, and the components of 
f may be deduced from equations (3)- (5). Letting tildes denote nondimensional quantities 
(e.g. R = R/R*), and using GM, = m a 3 t, we obtain 

u K = m\l*(V) X, \ (62) 

Similarly, the quantity ( is given by 

r 3 

( = ml ^. (63) 

16 R 

When the centrifugal radius equals the stellar radius, £ = 1 at the surface of the star. From 
the above equation, the corresponding nondimensional time r is 2.54. The infall and wind 
densities in nondimensional form are 

Pi = J^, (64) 
R z uk 

a 

Pw = 5^-, (65) 
R z v 

where the infall density shape function is 

S' 1 EE -f r [1 + 2CP 2 (COS0 O )]. (66) 

Dropping the tilde notation, the dimensionless form of the dynamical equations is 

~- u, (67) 



Dt 

D(\n(Aa)) 1 
Dt ~ R 2 a 



a(cos 7 -^)-^(/ n -^) 
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Du 



(68) 



Dt R 2 a 



«(co S7 -^)<-^( /n -^)u: 



~m ^ 2 . (69) 

Note again that only three parameters, the dimensionless ratios a, r, and z/, enter the final 
equations. 

In order to launch a shell, the driving wind speed must be sufficient to yield a shell 
velocity u r > 0, which in turn implies that the wind ram pressure exceed the infall ram 
pressure at the stellar surface. In terms of our chosen units, we require v > u ram , where 

v ram — 71 : ^ TToT (Jfy 

a (1 + m% T d /8) 

is the value of v that balances ram pressures at the stellar pole, and v esc = v / 2m r 
(c.f. eq.[62]) is the escape speed at the stellar surface. 



6. Numerical Results 

6.1. Breakout versus Recollapse 

Before discussing the details of our numerical results, we note that the mass of the shell 
typically will be dominated by the swept-up infall. In the limit of a stationary, spherical 
shell, the ratio of wind to infall contributions is a (< 1), but outward motion of the shell 
further increases the fractional contribution of infall to the shell mass. Although our 
calculations conserve the z-angular momentum of the incident, infalling gas, the centrifugal 
support of the shell is modest. The primary effect of infall angular momentum is to establish 
an asymmetric flow, which is then swept up by the expanding shell. 

A sample evolutionary sequence is shown in Figures 5-7 for a typical choice of 
nondimensional parameters (a = 1/3, r = 4, v — 4.7). Here, the results are displayed 
both nondimensionally and in physical units. For the latter, we assumed standard values 
of the dimensional quantities: R* = 3i? , Q = 2 x 10~ 14 s _1 , a Q = 0.2 km s _1 . These 
choices imply an accretion rate of 1.85 x 10~ 6 M yr _1 , a time since protostar formation of 
38,000 years, and a wind speed of 159 km s _1 . The shell first elongates as it expands up 
to and beyond the disk radius, indicated in Figure 5. This outward motion never ceases, 
as the wind is able to drive back the infalling envelope in all directions. We refer to such 
an outcome as breakout. While at intermediate times (see Figure 6) the shell is quite 
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elongated and a bipolar geometry is obtained, at late times (Figure 7) the shell takes on a 
more spherical appearance. In this spatial regime, well beyond the disk radius, the ambient 
density is also spherically symmetric. 

Figure 8 shows that dramatically different results are obtained with very similar input 
parameters (a = 1/3, r = 4, v — 4.4). Dimensionally, the wind speed has now been lowered 
to 148 km s _1 . The shell initially rises as before, but it stalls and falls back to the star 
at nearly free-fall speed. As the shell plummets to the stellar surface, initial ripples are 
amplified; higher numerical resolution would be needed to follow this apparent instability. 
We will not be interested in such details, but only in whether or not the shell succeeds in 
driving back the infall. Shells that break out are typically quite smooth and regular. Even 
for those that recollapse, the structure is smooth up to the point of turnaround. Note that 
recollapse is caused primarily by gravity acting on the heavy shell, rather than infall ram 
pressure. Turning off gravity in the code, breakout is achieved for a wind speed less than 
half of that in Figure 8. 

For those shells that do achieve breakout and have traveled far relative to the initial 
stellar radius, all quantities follow a power-law determined by the slope in the infall density 
law. For R cen <C r <C a t, pi oc r~ 3 / 2 , so the swept-up mass varies as r 3 / 2 , while the 
shell's momentum grows linearly with time. The resulting momentum-driven snowplow has 
R oc t 4 / 5 and V r oc r 1 / 5 . 

6.2. Critical Solutions 

We have mapped the parameter space of solutions as a function of a, r, and v. For a 
given value of a, the r — v plane is divided into three regions: one in which shells cannot 
even advance beyond the stellar surface, one in which the shell may initially rise, but is 
pulled back, and one in which the shell breaks out. These solution regimes will be called the 
crushed wind, the trapped wind, and the escaping wind, respectively. In the r — v plane, 
the boundary between the crushed wind and the trapped wind is the locus u ram (a, r) given 
by equation (70), where infall and wind ram pressures balance at the stellar surface. We 
similarly denote by u crit the minimum nondimensional wind speed necessary for breakout. 
Figure 9 plots this quantity as a function of time, for three representative values of a. 
Again, we display the results both nondimensionally and in physical units, employing our 
fiducial input parameters. Table 1 also lists values of v crit at selected times. Returning 
to the figure, we see that, at very early times, v cr u increases as r 1 / 2 , and is also inversely 
proportional to a. At intermediate times such that R cen ~ R*, the centrifugal deflection 
of infalling gas lowers the infall density, and u crit consequently falls. Once again, u crit at 
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any time during this epoch is proportional to a -1 . Throughout early and intermediate 
times, whether or not the shell breaks out is determined by the product a u, that is, on the 
momentum loss rate of the wind. Finally, v crit again rises at late times. In this case, the 
critical wind speed is independent of a and increases as t 1 / 2 . 

Let us see in more detail how these results arise. The wind must combat the infall 
ram pressure and the gravitational force on the shell. At early times, we may neglect 
centrifugal distortion of the infall and consider only the spherically symmetric problem. In 
this case, the infall ram pressure at fixed radius varies as r 1 / 2 . (Recall equation (39).) The 
gravitational force per unit shell mass rises as r 1 . To determine the force per area, we first 
note that the infall density pi falls as r -1 / 2 , a consequence of the rising freefall velocity. For 
a light wind, p w <^ Pi, we may neglect the mass contribution of the wind. Since most of the 
mass of the shell is swept-up, a oc pi oc r -1 / 2 , (c.f., equations (62), (64)) so the gravitational 
force per unit area on the shell aV also behaves as t 1 / 2 . Comparing the wind ram pressure 
oc a v to the sum of infall ram pressure and gravity, we obtain v crit oc r 1 / 2 /a. 

The downturn in u crit at r ~ 1 occurs because the star is now interior to the centrifugal 
radius. The decrease in the infall density along the z-axis relative to spherical accretion 
implies that breakout becomes easier. Because the shell's mass is dominated by the infall 
contribution, the decrease in the infall density along the axis decreases both the infall ram 
pressure and the gravitational force on the shell. At late times, the centrifugal radius has 
grown so large that the primary source of mass input to the shell is the wind. In this 
regime, the wind speed necessary for breakout is proportional to the stellar escape velocity, 
which in turn varies as r 1 / 2 . 

We compare these numerical results for v crit with three analytical approximations in 
Figure 10. A hypothetical wind of speed v = v esc (the diagonal, dash-dot line) would be 
fast enough to break out if we neglect the mass and momentum fluxes to the shell from 
infall. We see that except for late times (r > 10), the wind must be considerably faster 
than the escape speed to drive a swept-up shell that will break out. At late times, our 
constant speed wind can break out with v < u esc because we have neglected gravitational 
deceleration of the preshock wind. 2 The second approximation is that of u = u rami shown 



2 Using modified wind conditions of such a coasting, decelerating wind of constant specific energy 
e = — GM<,/r, and numerically determining the corresponding critical speed for breakout with our 

code (filled triangles in Figure 10), we see that the late-time critical curves are shifted vertically and converge 
to the escape speed. At early times, since breakout speeds were already much larger than the escape speed, 
the difference between a constant speed and a decelerating wind is negligible, and the critical curves are 
unchanged. Because no observed winds actually decelerate, we consider the constant velocity curves to be 
more applicable. 
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by the dashed line, yielding ram pressure balance of wind and infall at the stellar surface. 
This approximation underestimates the necessary wind speed for breakout by a factor of 1.9 
at early times, when the infall is spherically symmetric, but it underestimates the critical 
speed by a much larger factor once the infall becomes significantly aspherical. 

The two foregoing analytical approximations give a poor agreement with v crit because 
the first neglects the infall ram pressure, while the second neglects the gravitational force on 
the shell. To include in an analytic fashion the effect of both infall ram pressure and gravity 
on the shocked gas within the shell, we specify the condition that the initial postshock gas 
be rising at the escape speed. According to equation (50), in order to have u , r = V esc , we 
obtain the necessary condition 

v = v esc {l + 2n 1 ' 2 ). (71) 

We cannot immediately apply this equation, because n = pi/p w itself depends on the 
wind speed. The density ratio is found from equations (64)-(66), where for 9 = we have 
f r = —\/2 and P 2 (cos9 ) = 1, which yields 

V = b 2 — , (72) 

L'esc 

where we have defined for brevity 

b = 1/^(1 + 2 Co). (73) 

The solution v' to equations (71), (72) is found as the root of a quadratic equation, giving 

u' = fu esc [b+Vb^Tl] 2 , (74) 

which for / = 1 is the value of v required to give the postshock, mixed gas escape speed at 
the stellar surface. The ad hoc factor / allows us to consider the wind to be some fraction 
/ of the speed required to give the shell escape speed at launch. The corresponding curve 
for a = 0.1 and / = 1 is shown in Figure 10 (dotted curve). If we neglect changes in the 
infall momentum flux per unit solid angle once the shell is launched, v 1 is expected to be an 
estimate of the necessary wind speed to launch a shell that will break out. In practice, as the 
shell advances, the forces on the shell change, and typically the wind weakens less quickly 
than the infall, so the above condition overestimates the required wind speed. Except for 
the vertical offset, however, the curve follows closely the shape of the numerically-derived 
critical curve, suggesting / ps constant throughout the evolution. By choosing / = 0.45, 
so as to fit the late-time behavior, an approximate fit to the critical curves is shown in 
Figure 9. 
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6.3. Generalization to Anisotropic Wind 

The previously described calculations may be extended to anisotropic winds by denning 
angular dependences of the wind mass and momentum fluxes according to 

p w V w = -^f w (0), (75) 
^ = ^T9M, (76) 

where V w is the streamline-averaged wind speed and the functions f w and g w are normalized 
to have unit average value over An steradians. Although a diversity of shell shapes 
may be generated in this fashon, we focus only on the behavior at the symmetry axis, 
where breakout is easiest. If the properties of the wind are smooth near the pole, i.e. 
/t«(0) = fltu(O) = 0; the shell has R'(0) = and looks locally as though it is driven by an 
isotropic wind. 

Figure 6 shows a numerical example. Here the dashed curves represent a shell driven 
by an anisotropic wind specifically chosen to have the same mass and momentum fluxes 
towards 9 = as the isotropic wind (solid curves). We have chosen the asymmetric driving 
wind to have a — 1/6 and a density p w oc f w — g w — |(1 + 3 cos 2 9), with V w independent 
of 9, while the values of v and r are unchanged from the spherical wind case. This angular 
dependence represents the lowest-order expansion that is axisymmetric and of even parity, 
and has been chosen arbitrarily to give a pole-to-equator density contrast of 4. Because 
the derived polar behavior R (t) of a shell driven by an anisotropic wind follows that of a 
spherical wind having the same polar mass and momentum fluxes, the critical wind speed 
is crit for breakout of an anisotropic wind may be determined from our numerical curves. 
Defining v = V w /V*, the equivalent spherical wind corresponds to a sp h = ex f w (0) and 
v sp h = v g w (0)/ f w (0). A wind that is focused preferentially towards the poles mimics a 
spherical wind corresponding to a higher value of a and a reduced value of u critl breaking 
out at an earlier time (see below). In our example, the a — 1/6 wind has a sp h = 1/3 and 
the shell propagates along the axis like a bubble from an isotropic wind of higher a-value. 



6.4. The Breakout Time 

Given the curves for the critical wind speed for breakout (Figure 9), we now wish to 
estimate the time at which breakout occurs. We first note that a wind whose speed is 
independent of r would experience no trapped phase, since the curves for u crit decrease as 
t — > 0. Although we are not offering a general theoretical account of protostellar winds, we 
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note that a constant wind speed is unlikely on purely empirical grounds. Stars of all masses 
and ages tend to have V w -vahies roughly equal to the appropriate surface escape speed 
(Lamers & Cassinelli 1999). We therefore recast Figure 9 as a plot of V w /V esc = vjv esc 
versus r (Figure 11). For comparison, we have also plotted v ram jv esc ^ to delimit the end 
of the crushed wind phase, using equation (70) for v ram . Suppose, then, that vjv esc were 
strictly a constant, equal to unity. Taking the case a — 1/3 in Figure 11, we see that 
the wind would first be able to advance from the stellar surface at t ~ 25, 000 years, but 
breakout would be achieved at t ~ 54,000 years. For the a = 1/10 model, these times 
increase to 41, 000 and 98, 000 years, respectively. Our numbers rely on a relatively modest 
value of f2 Q = 2 x 10 _14 s _1 . Based upon equation (52), increasing f2 G by a factor of two 
would decrease the above breakout times by a factor of 0.63. Thus, the trapped outflow 
phase may be a significant fraction (~ 60%) of the time prior to breakout. 

For the example we have shown, when the shell breaks out, it does so at essentially 
all angles. However, this is not the case for all runs and is not true in general. Thus, our 
breakout time is to be interpreted as the time at which the wind first breaks out along the 
z— axis. At larger angles, the wind may be unable to break out until a later evolutionary 
time, resulting in simultaneous outflow and infall. Usually, however, the wind escapes 
globally, as shown in Figures 5-7. The shell aspect ratio then begins and ends as unity, 
either at the stellar surface or far from the star. The maximum distortion, which is still 
modest, is attained at intermediate times, when the shell is a few AU from the stellar 
surface. 



7. Discussion 

7.1. Behavior during the Trapped Wind Phase 

Our recollapsing shells generally become unstable during the trapped wind phase. 
Portions that are initially inside neighboring regions fall faster, increasing any perturbation 
of an initially smooth surface. While our calculations do find this strong dynamical effect, 
we cannot detect the nonlinear instability of Vishniac (1994), which relies on a finite shell 
thickness. Our results make it clear, in any case, that the true behavior during the trapped 
wind phase is exceedingly complex. After a portion of the shell collapses to the stellar and 
disk surfaces, a new shock is immediately driven outward by the steady wind. Thus, there 
will simultaneously be contracting and expanding patches at different angles. Portions of 
the wind may escape temporarily between recollapsing fragments before falling back. While 
the details of this process cannot be followed by our computational method, we believe 
our calculated breakout times should still be reasonably accurate. In any event, this time 
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cannot decrease by more than about a factor of two, since the infall crushes the wind 
entirely at an earlier epoch. 

One observational signature of the trapped wind phase would be fluctuations in 
luminosity due to the rising and recollapse of the shell. This luminosity arises from 
both the shocking of the wind and ambient gas. When a portion of the shell is rising, 
the corresponding wind component of luminosity decreases because of the lowered shock 
velocity. The infall component of the luminosity may either increase (because the shock 
velocity is raised), or decrease (because |uj| = u esc decreases at larger radius.) All luminosity 
fluctuations should increase in both period and amplitude with evolutionary time, because 
as breakout is approached, the shell rises higher before recollapse. Further study is needed 
to determine the observed magnitude of these fluctuations, as well as the transport of 
radiation through the infalling gas. As indicated earlier, the wind may also break out 
along the axis, while portions of the shell at larger angles continue oscillatory behavior. In 
this case, one would observe a bipolar jet accompanied by quasi-periodic fluctuations in 
luminosity. 



7.2. Comments on Collimation 

Figures 5-7 show that the lobe-like appearance of our shells is a transient phenomenon 
created by anisotropy in the infalling gas. In particular, this elongated morphology cannot 
be identified with the collimation seen in Herbig-Haro jets. Our assumed mixing of shocked 
wind and infall is, in fact, an agent for decollimation. A fluid element of the wind, initially 
moving radially outward, acquires a positive ^-velocity. Thus, there is a tendency for gas to 
move towards the disk plane, despite the (temporary) elongation of the shell. 

Our mixing assumption was made purely for computational simplicity. In reality, there 
must be considerable shear between the two shocks. This shear may have several important 
consequences. First, shocked wind may be refracted along the shell's inner surface toward 
the polar axis, resulting in a more jet-like flow. Such a geometry would resemble the early 
outflow model of Canto (1980), although the latter pictured a shell in pressure balance 
with a static, external medium. A second effect is that shocked, ambient gas will move 
tangentially toward the equator. The overall result will be a shell of lower mass that is 
less confined by stellar gravity. The duration of the trapped wind phase will therefore be 
reduced, resulting in earlier breakout. However, we note that even in the limit of no mixing, 
the wind must still be crushed by infall at a sufficiently early epoch. In future calculations, 
we hope to explore quantitatively the consequences of relaxing the mixing assumption. 
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In a series of numerical simluations carried out through the 1990s, A. Frank and 
colleagues found that the interaction of a spherical wind and infall produces very strong 
collimation. However, these studies (beginning with Frank & Noriega-Crespo 1994) modified 
the equations for infall, thus obtaining a much larger density anisotropy than adopted here. 
For example, Frank & Mellema (1996) used an equator-to-pole density contrast of 50 — 70, 
while Delamarter et al. (2000) allow this ratio to exceed 10 3 . For comparison, the density 
contrast in our outflow only exceeds 50 only for 0.96 £ R/R cen £ 1.08. Our escaping shell 
spends relatively little time within this region, sweeping up a small amount of mass. We 
conclude that a spherical or modestly anisotropic wind undergoes little collimation through 
interaction with a physically realistic infall, unless it is through the shear effect described 
above. 

Our finding of an early, trapped wind phase is a similarly robust result, although it 
was missed in previous simulations. Frank & Mellema adopt such a large wind velocity 
(V w > 500 km s" 1 ) that the infall ram pressure is overwhelmed. Mellema & Frank (1997) 
did find oscillating shells that, under some conditions, collapsed. However, the driver of 
their oscillations was a variable wind, so that collapse always coincided with the weakest 
outflow phase. More recently, Delamarter et al. (2000) utilized as an infalling background 
the flattened, rotating cloud of Hartmann, Calvet, & Boss (1996). Although they claim 
that especially weak winds may be stifled entirely by infall, their numerical results always 
show the wind escaping through at least a narrow solid angle through the symmetry axis. 

We maintain, based on the arguments presented above, that early wind crushing and 
eventual breakout are inescapable features within a realistic account of the wind-infall 
interaction. In contrast, a strong, untrapped wind that is present ab initio, as in previous 
simulations, would have already modified substantially the background infall, rendering 
subsequent results of dubious validity. Our own calculations of breakout can provide useful 
starting conditions for future investigators who wish to follow the wind evolution well 
beyond this critical, early phase. 

7.3. Further Observational Considerations 

We have seen that, for canonical values of our parameters, the wind is trapped for 
a period of roughly 50, 000 yr. Although the evolutionary status of the most embedded 
sources is by no means secure, those designated as Class (Andre et al. 1993) are generally 
considered to be the youngest. They are detectable only at far-infrared and submillimeter 
wavelengths, and have spectral energy distributions corresponding to dust temperatures 
of roughly 30 K. An additional argument for their extreme youth is their relatively low 
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population. In a submillimeter survey of p Ophiuchi, Motte et al. (1998) found that Class I 
sources outnumber those of Class by a factor of 10. Since the former have traditionally 
been assigned ages of order 10 5 yr (Kenyon et al. 1990), Class stars could be as young as 
10 4 yr, if the population is forming in steady-state fashion. Visser et al. (2002), in their 
recent SCUBA survey of Lynds dark clouds, have challenged these figures. They find nearly 
equal numbers of Class and I sources, and conclude that the former also have ages of 
roughly 10 5 yr. 

Regardless of the outcome of this controversy, both surveys, as well as other studies, 
have found that essentially all stars with massive, dusty envelopes produce winds, as 
evidenced by associated CO outflows. Where, then, is the trapped phase? We see two 
possible answers to this question. The first, and less likely, is that our choice of input 
parameters requires adjustment. However, the necessary changes would be severe. For 
fixed a, v, and r, tb re ak scales as rV 3 a" 1 / 3 Q _2//3 . It strains credibility to suppose, for 
example, that the cloud sonic speed a is so high that U rea k is reduced by an order of 
magnitude. In addition to this change of scale, we may change nondimensional solutions, 
by assuming a faster wind (greater v). As seen in Figure 11, the breakout time is relatively 
insensitive to v until a value is reached where the trapped wind disappears entirely (where 
the dashed curves are horizontal). This option, however, would imply very high wind 
speeds (V w ~ 3V esc for a = 1/3, or V w ~ 10 V esc for a = 1/10). As noted in §6.3, making 
the wind more anisotropic effectively raises a and also leads to earlier breakout. Winds 
launched via the magnetocentrifugal mechanism, however, typically have modest anisotropy 
until they propagate a considerable distance (see, e.g., Najita & Shu 1994). A second, and 
more plausible, answer is that we are already witnessing at least partially trapped winds. 
Following the arguments in §7.1, it may be that simultaneous infall and outflow occur 
well before our idealized model indicates the onset of breakout. It would be extremely 
interesting, in this regard, to monitor the temporal variability of outflows from Class 
sources. 



8. Conclusions 

We have presented a detailed numerical calculation of the interaction between a 
spherical, protostellar wind and an anisotropic, infalling envelope. The anisotropy in our 
model is derived from the rotation of the star's parent cloud core. Collapse is assumed 
to proceed in inside-out fashon from a singular, isothermal sphere, in the absence of 
magnetic forces. We have idealized the double-shock interaction region as a thin shell that 
is well-mixed internally. After demonstrating that our previous, quasi-steady solutions 
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are dynamically unstable, we have followed the shell dynamics in a fully time-dependent 
manner. 

At very early times, the wind is crushed by ram pressure associated with infalling 
matter. Wind breakout is delayed by the gravitational force exerted on the shocked 
material. The shell must not only be supported by wind ram pressure, but must have 
sufficent kinetic energy to escape the star's potential well, even while gathering additional 
mass from the envelope. Breakout would occur earlier if either the wind or infall were more 
anisotropic, or if we were to include the internal shear between the shock surfaces. 

Bearing in mind the observations of jets and outflows from embedded stars, we 
acknowledge two important limitations of our model. First, our assumption of complete, 
early trapping may be incorrect, because of shell fragmentation. Second, the observed 
collimation of stellar jets over long distances is not approached asymptotically in our 
calculation. Instead, our shells, while they are temporarily elongated, eventually become 
spherical. Jet collimation may be a consequence of magnetic pinching in the wind itself, or 
of a cloud background that is very different from the one we assume here. An initially more 
anisotropic cloud core will, of course, yield an altered pattern of infall. Alternatively, a 
shallower falloff in the background density will give rise to crossing shocks in the interaction 
region (see, e.g., Canto, Raga, & Binette 1989). These shocks further aid in wind 
collimation. 

FPW acknowledges helpful discussion with S.Lizano, and is grateful to the NSF 
International Researchers Fellowship Program and CONACyT/Mexico for financial 
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Table 1. Critical Wind Speed v crit 
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Fig. 1. — Element of the shell, with normal and tangential directions shown. The angle 
7 is measured clockwise from the radial to the normal direction. The shell thickness is 
An, while the arc length at constant azimuthal angle is As. The azimuthal width of the 
patch is Aw = RsmQAcp. The length of the chord slicing the shell at constant (0,4>) is 
Ar = An sec 7. 

Fig. 2. — Approach of the shell to the steady-state bow shock solution. In the frame of the 
star, located at the origin, the ambient medium moves with velocity —V*z. The initially 
spherical shell began at 3 R sun , much smaller than the standoff radius R = 2.45 x 10 17 cm. 
Solutions are shown for elapsed time increasing by a factor of two from ^ to 4 times the 
crossing time R /V*. Triangles denote gridpoints, while the analytic, steady-state solution 
is given by the solid curve. 

Fig. 3. — Fractional errors for the bow shock test case after reaching near-equilibrium. 
Because the normal velocity vanishes in steady state, the quantity u n /u t refers to the normal 
velocity of the gridpoints divided by the analytic solution for the tangential velocity. 

Fig. 4. — Instability of an Initial Steady-State Solution: Polar radius vs. time for shells 
starting close to the steady-state solution of Paper I. Shells begun at slightly smaller radius 
collapse, while shells begin at slightly larger radius expand indefinitely. The equilibrium 
radius is shown (boldface) for the solution corresponding to a — 0.1, ( = 2.45. 

Fig. 5. — Time-evolution of an escaping shell (early evolution), corresponding to a — 1/3, 
v = 4.7, and r = 4.0. The protostellar age since core formation is 3.8 x 10 4 years. The 
shapes correspond to equal time intervals of 0.016 years. The subsequent evolution of this 
shell is shown in the next two figures on a larger scale. The scale of the centrifugal radius 
is indicated by the disk. Lengths are in cm on the left and bottom axes, and in units of the 
stellar radius on the top and right axes. 

Fig. 6. — Time-evolution of an escaping shell (further evolution), corresponding to the same 
parameters as Figure 5 but for with increasing times steps and on a larger scale. The 
innermost curve of this figure is the same as the outermost of that Figure. However, the 
elapsed time increases by a factor of two with each curve: 4,8,16,32, and 64 time units. The 
dashed curves represent an "equivalent" model driven by an asymmetric wind (See §6.3) 

Fig. 7. — Time-evolution of an escaping shell (late time), corresponding to the same 
parameters as Figures 5 and 6. Here the shell has gone well beyond the centrifugal radius 
of the infall and is in the increasingly spherically symmetric infall region. Hence the shell 
becomes more spherical as it sweeps up this material. The elapsed time doubles with each 
successive curve. The innermost here corresponds to 128 time units of Figure 5, while the 
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outermost is greater by a factor of 64 (8192 time units), which corresponds to 134 years since 
launch. 

Fig. 8. — Time-evolution of a recollapsing shell corresponding to a — 1/3, v — 4.4, and 
r = 4.0. The critical wind speed for these parameters is u crit = 4.48. Solid curves show 
the rising phase, while dotted curves display the subsequent recollapse). Shown for equal 
time intervals of 1/4 the time to reach the highest point along the z-axis, which is the same 
interval as the unit used in Fig. 5. 

Fig. 9. — Minimum breakout wind speed versus evolutionary time. The three loci (solid 
curves) correspond to three ratios, a, of the wind mass loss to infall accretion rate. For a 
given a, the region above the curve corresponds to breakout, while that below the curve 
corresponds to recollapse. An analytic fit to the numerical solutions is shown in dashed 
curves. 

Fig. 10. — Minimum breakout wind speed versus evolutionary time, compared with analytical 
arguments. All curves correspond to a — 0.1. At bottom left (dashed curve), the ram 
pressure balance condition at the launch point. Bottom right (diagonal, dot-dashed line), 
the condition of initial wind speed equal to the escape speed. Curve labeled v crit : numerical 
results for the critical wind speed. Top curve (dotted): the condition that the shocked shell 
initially move at the escape speed. Filled triangles: critical condition for a decelerating wind 
(see text). 

Fig. 11. — Critical wind speed for breakout (solid curves), in units of the free-fall (escape) 
speed, as a function of evolutionary time. The corresponding a-values are shown, as well 
as the wind speed necessary for ram pressure balance at the stellar surface (dashed curves). 
Assuming wind launch conditions v w jv e sc = 1 (i-e. following the horizontal line in this 
figure), evolution begins at the left edge of the plot with the wind unable to advance beyond 
the stellar surface until the line intersects the appropriate dashed curve. Then the trapped 
wind phase lasts until the line intersects the corresponding solid curve for breakout. 
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